4  CVM ppp

4.1 Pretratamientos de la Información

Cargar Recursos Externos

source("scripts/recursos.R")
source("scripts/funciones.R")

Definición de variables generales

## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <- 
  "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Definición de variables específicas

# List of crimes/delitos and cities/ciudades
crimes = c(
  "Minor property offenses",
  "Robberies",
  "Burglaries",
  "Drunkennes, damages and disorders",
  "Family violence and aggression",
  "Injuries, drugs and weapons",
  "High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")

#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]


Lectura de Datos Espaciales

Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.

city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))

El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame.


Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)

# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)


Ventana de Trabajo para ppp (Point Pattern Plane)

Crea un objeto de clase “ppp” que representa un conjunto de datos de patrones de puntos en el plano bidimensional. Por lo anterior, se debe crear una ventana de observación w.

#make windows
  ext <- raster::extent(cityp)
  x_min <- ext[1] + 100
  x_max <- ext[2] - 100
  y_min <- ext[3] + 100
  y_max <- ext[4] - 100
  w <- as.owin(c(x_min,x_max, y_min, y_max))


Transformaciones a los datos.

Se utilizó funciónes de transformacion general mencionadas en Section 6.2 que corresponden a las siguientes:

  • Transformación Cúbica
  • Transformación Logarítmica
  • Cáculo del Lag Espacial
cityp <- cityp %>%
    st_as_sf(crs = 4326) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_cub, .names = "cub{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_log1p, .names = "log{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = lag_spatial, nb, .names = "lag{.col}")) %>%
    as("Spatial")
zona id x y pflot viol crim abus amen robo hurto asalto denspob educ tamhog constr ofcom vulzon mobil lpflot neduc ldenspob lofcom lconstr cubhurto cubrobo cubasalto cubamen cubabus cubcrim cubviol loghurto logrobo logasalto logamen logabus logcrim logviol laghurto lagrobo lagasalto lagamen lagabus lagcrim lagviol
4101161005 22 -71.239 -29.862 1.000 0 0 0 0 0 0 0 144.756 10.867 3.170 4328.618 0.000 0.367 0.079 0.178 0.477 0.767 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.417 0.250 0.083 0
4101161005 32 -71.243 -29.863 1.000 0 0 0 0 0 0 0 227.772 10.716 3.393 3830.009 0.000 0.367 0.079 0.178 0.466 0.846 0.000 0.829 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.083 0.917 0.583 0.083 0
4101161005 35 -71.240 -29.863 1.034 0 0 0 0 0 0 1 171.133 11.801 3.172 4328.618 0.000 0.367 0.079 0.182 0.550 0.796 0.000 0.840 0 0 1 0 0 0 0 0 0 0.693 0.000 0.000 0 0 0.000 0.083 0.000 0.500 0.250 0.083 0
4101161005 36 -71.239 -29.863 1.000 0 0 0 0 0 0 0 154.603 10.870 3.135 4328.618 0.000 0.367 0.079 0.178 0.478 0.779 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.417 0.250 0.083 0
4101161005 37 -71.238 -29.863 1.025 0 0 0 0 0 0 0 155.145 10.858 3.148 3230.309 0.000 0.367 0.079 0.181 0.477 0.779 0.000 0.815 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.083 0.083 0.083 0.000 0.250 0.083 0
4101161001 48 -71.244 -29.864 1.125 0 0 0 0 0 0 0 168.179 11.958 3.336 2384.222 14.682 0.292 0.122 0.193 0.562 0.793 0.574 0.788 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.000 0.250 0.500 0.417 0.083 0
4101161005 49 -71.243 -29.864 1.071 0 0 0 0 0 0 0 212.743 10.978 3.415 3795.759 0.000 0.367 0.079 0.187 0.486 0.834 0.000 0.828 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.250 0.917 0.583 0.083 0
4101161005 50 -71.242 -29.864 1.000 0 0 1 1 0 0 0 220.188 10.705 3.409 4026.613 0.000 0.367 0.079 0.178 0.465 0.840 0.000 0.834 0 0 0 1 1 0 0 0 0 0.000 0.693 0.693 0 0 0.000 0.083 0.000 1.000 0.417 0.083 0
4101161005 51 -71.241 -29.864 1.021 0 0 1 1 0 0 0 243.501 10.120 3.337 4328.618 0.000 0.367 0.079 0.180 0.419 0.857 0.000 0.840 0 0 0 1 1 0 0 0 0 0.000 0.693 0.693 0 0 0.000 0.083 0.083 0.667 0.250 0.000 0
4101161005 52 -71.240 -29.864 1.022 0 0 0 0 0 0 0 320.023 10.720 3.324 4328.618 0.000 0.367 0.079 0.180 0.466 0.905 0.000 0.840 0 0 0 0 0 0 0 0 0 0.000 0.000 0.000 0 0 0.000 0.083 0.167 0.583 0.250 0.000 0


4.2 Evaluación Cramer Von Mises

4.2.1 Modelos MLB

Para transformar los datos espaciales de los delitos y las predicciones de los modelos se utilizó la función propia ppp_maker() Section 6.4 que retorna el objeto ppp.

Lectura del Modelo.

model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))

Transformaciones de las predicciones del modelo a ppp.

ppp_mbl_pred <- ppp_maker(model = model_mlb, city = cityp, 
                     delito = delito, ppp_predict = T, w = w)

plot(density(ppp_mbl_pred, adjust=.25))

Transformaciones de los delitos reales a a ppp.

ppp_mbl <- ppp_maker(model = model_mlb, city = cityp, 
                          delito = delito, ppp_predict = F, w= w)

plot(density(ppp_mbl, adjust=.25))

Evaluación cdf.test.()

cvm_mlb <- cdf.test(ppp_mbl, as.im(ppp_mbl), test="cvm")
plot(cvm_mlb)

attr_model <- attr(cvm_mlb, which = "frame")
im <- attr_model$values$Zimage 
r_mlb <- raster(im)
proj4string(r_mlb)=crs_utm
values(r_mlb) <- ifelse(values(r_mlb)==0, NA, values(r_mlb))
plot(r_mlb)

pred_mlb <- predict_df(model = model_mlb, city = cityp, longlat = F)
pred_mlb <- st_as_sf(x = pred_mlb,                         
                       coords = c("x", "y"),
                       crs = 32719) %>% 
  filter(!is.na(delitos))

mapview(r_mlb, na.color= NA)+
  mapview(pred_mlb, zcol = "ypred",hide = TRUE,  cex = 2)+
  mapview(pred_mlb, zcol = "delitos", hide = TRUE,  cex = 2)
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
+wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

4.2.2 Modelos SAR

city <- ciudades[obs] 
model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))

ppp_sar_pred <- ppp_maker(model = model_sar, city = cityp, 
                            delito = delito, ppp_predict = T,
                            w = w, nb = nb)
This method assumes the response is known - see manual page
plot(density(ppp_sar_pred, adjust=.1))

ppp_sar <- ppp_maker(model = model_sar, city = cityp, 
                       delito = delito, ppp_predict = F, 
                       w = w, nb = nb)
  
plot(density(ppp_sar, adjust=.25))

# st <- syrjala.test(ppp_mbl, ppp_mbl_pred, nsim = 999)
cvm_sar <- cdf.test(ppp_sar, as.im(ppp_sar_pred), test="cvm")
plot(cvm_sar)

attr_model <- attr(cvm_sar, which = "frame")
im <- attr_model$values$Zimage 
r_sar <- raster(im)
proj4string(r_sar)=crs_utm
values(r_sar) <- ifelse(values(r_sar)==0, NA, values(r_sar))
# plot(r_sar)
pred_sar <- predict_df(model = model_sar, city = cityp, longlat = F)
This method assumes the response is known - see manual page
pred_sar <- st_as_sf(x = pred_sar,                         
                       coords = c("x", "y"),
                       crs = 32719) %>% 
  filter(!is.na(delitos))

mapview(r_sar, na.color= NA)+
  mapview(pred_sar, zcol = "ypred",hide = TRUE,  cex = 2)+
  mapview(pred_sar, zcol = "delitos", hide = TRUE,  cex = 2)